//----------------------------------------------------------------------------
//  XC program; finite element analysis code
//  for structural analysis and design.
//
//  Copyright (C)  Luis Claudio Pérez Tato
//
//  This program derives from OpenSees <http://opensees.berkeley.edu>
//  developed by the  «Pacific earthquake engineering research center».
//
//  Except for the restrictions that may arise from the copyright
//  of the original program (see copyright_opensees.txt)
//  XC is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or 
//  (at your option) any later version.
//
//  This software is distributed in the hope that it will be useful, but 
//  WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details. 
//
//
// You should have received a copy of the GNU General Public License 
// along with this program.
// If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------------
/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 2001, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** Reliability module developed by:                                   **
**   Terje Haukaas (haukaas@ce.berkeley.edu)                          **
**   Armen Der Kiureghian (adk@ce.berkeley.edu)                       **
**                                                                    **
** ****************************************************************** */
                                                                        
// $Revision: 1.3 $
// $Date: 2004/01/22 01:56:47 $
// $Source: /usr/local/cvs/OpenSees/SRC/reliability/analysis/transformation/NatafProbabilityTransformation.cpp,v $


//
// Written by Terje Haukaas (haukaas@ce.berkeley.edu)
//

#include <reliability/analysis/transformation/ProbabilityTransformation.h>
#include <reliability/analysis/transformation/NatafProbabilityTransformation.h>
#include <reliability/domain/components/RandomVariable.h>
#include <reliability/domain/components/CorrelationCoefficient.h>
#include <reliability/domain/distributions/NormalRV.h>
#include <utility/matrix/Vector.h>
#include <utility/matrix/Matrix.h>
#include <reliability/analysis/misc/MatrixOperations.h>
#include <cmath>
#include <cstring>


XC::NatafProbabilityTransformation::NatafProbabilityTransformation(ReliabilityDomain *passedReliabilityDomain,
                                                                                         int passedPrintFlag)
:ProbabilityTransformation()
{
        theReliabilityDomain= passedReliabilityDomain;
        printFlag= passedPrintFlag;


        // Find and set problem size (number of random variables)
        nrv= theReliabilityDomain->getNumberOfRandomVariables();


        // Create/initialize vectors and matrices
        x= new Vector(nrv);
        u= new Vector(nrv);
        jacobian_x_u= new Matrix(nrv,nrv);
        jacobian_u_x= new Matrix(nrv,nrv);
        lowerCholesky= new Matrix(nrv,nrv);
        inverseLowerCholesky= new Matrix(nrv,nrv);
        correlationMatrix= new Matrix(nrv,nrv);


        // Establish correlation matrix according to the Nataf assumption
        setCorrelationMatrix(0, 0, 0.0);


        // Create object to do matrix operations on the correlation matrix
        theMatrixOperations= 0;
        theMatrixOperations= new MatrixOperations(*correlationMatrix);
        if (theMatrixOperations == 0) {
                std::cerr << "XC::NatafProbabilityTransformation::NatafProbabilityTransformation() - could " << std::endl
                        << " not create the object to perform matrix operations." << std::endl;
        }


        // Cholesky decomposition of correlation matrix
        int result= theMatrixOperations->computeCholeskyAndItsInverse();
        if (result < 0) {
                std::cerr << "XC::NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << std::endl
                        << " compute the Cholesky decomposition and its inverse " << std::endl
                        << " for the correlation matrix." << std::endl;
        }
        (*lowerCholesky)= theMatrixOperations->getLowerCholesky();
        (*inverseLowerCholesky)= theMatrixOperations->getInverseLowerCholesky();
}

XC::NatafProbabilityTransformation::~NatafProbabilityTransformation()
{
        if (correlationMatrix != 0) 
                delete correlationMatrix;
        if (lowerCholesky != 0) 
                delete lowerCholesky;
        if (inverseLowerCholesky != 0) 
                delete inverseLowerCholesky;
        if (jacobian_x_u != 0) 
                delete jacobian_x_u;
        if (jacobian_u_x != 0) 
                delete jacobian_u_x;
        if (x != 0) 
                delete x;
        if (u != 0) 
                delete u;
        if (theMatrixOperations != 0) 
                delete theMatrixOperations;
}



int 
XC::NatafProbabilityTransformation::set_x(Vector passedx)
{
        (*x)= passedx; // (later: check size of vector, etc.)
        return 0;
}

int 
XC::NatafProbabilityTransformation::set_u(Vector passedu)
{
        (*u)= passedu; // (later: check size of vector, etc.)
        return 0;
}

int 
XC::NatafProbabilityTransformation::transform_x_to_u()
{
        Vector z= x_to_z(*x);
        (*u)= (*inverseLowerCholesky) * z;
        
        return 0;
}


 XC::Vector 
XC::NatafProbabilityTransformation::meanSensitivityOf_x_to_u(Vector &x, int rvNumber)
{
        // Returns the sensitivity of 'u' with respect to [mean, stdv]
        // for the given random variable number (rvNumber)

        // Need four factors:

        // 1) (*inverseLowerCholesky)    DONE

        // 2) Vector z= x_to_z(x); 
        Vector z= x_to_z(x);

        // 3) DzDmean and DzDstdv= a vector of zeros and then:
        Vector DzDmean(x.Size());
        NormalRV aStandardNormalRV(1,0.0,1.0,0.0); 
        RandomVariable *theRV= theReliabilityDomain->getRandomVariablePtr(rvNumber);
        if(strcmp(theRV->getType(),"NORMAL")==0)
          {
            //double mu= theRV->getParameter1();
            double sigma= theRV->getParameter2();
            DzDmean(rvNumber-1)= -1.0 / sigma;
          }
        else if(strcmp(theRV->getType(),"LOGNORMAL")==0)
          {
            //const double lambda= theRV->getParameter1();
            //const double zeta= fabs(theRV->getParameter2()); // more here for negative lognormal?
            const double mean= fabs(theRV->getMean()); // more here for negative lognormal?
            const double stdv= theRV->getStdv();

            const double a= mean*mean+stdv*stdv;
            DzDmean(rvNumber-1)= 0.5*(-2.0*mean*mean*log(a)+4.0*mean*mean*log(mean)-3.0*stdv*stdv*log(a)+4.0*stdv*stdv*log(mean)+2.0*stdv*stdv*log(fabs(x(rvNumber-1))))/(pow(log(a)-2.0*log(mean),1.5)*mean*a);
                
//          double z= ( log ( fabs(x(rvNumber-1)) ) - lambda ) / zeta; // more here for negative lognormal?
//          double e= (stdv/mean)*(stdv/mean);
//          double d= 1 + e;
//          double f= 1.0 / (mean*d*zeta);
//          DzDmean(rvNumber-1)= f*(-d-e+e*z/zeta);
          }
        else if (strcmp(theRV->getType(),"UNIFORM")==0)
          {
                double pz= 0.39894228048*exp(-z(rvNumber-1)*z(rvNumber-1)/2.0);
                double a= theRV->getParameter1();
                double b= theRV->getParameter2();
                DzDmean(rvNumber-1)= -1.0/(pz*(b-a));
          }
        else
          {
            std::cerr << "WARNING: Cannot compute reliability sensitivity results for " << std::endl
                      << " type of random variable number " << rvNumber << std::endl;
            DzDmean(rvNumber-1)= 0.0; //aStandardNormalRV.getInverseCDFvalue(theRV->getCDFMeanSensitivity(x(rvNumber-1)));
          }

        // 4) The hardest part: DinverseLowerCholeskyDmean
        //                                                DinverseLowerCholeskyDstdv
        //    This is possible, but a bit tedious
        //    The procedure would be:
        //             R= L * L^T
        //            dR= dL* L^T + L * dL^T
        //            dR= dL* L^T + dL* L^T
        //            dR= 2 * dL* L^T
        //            dL= 0.5 * dR * L^(-T)
        //
        //    For now: do this by finite difference
        Matrix OrigInverseLowerCholesky= (*inverseLowerCholesky);
        RandomVariable *aRandomVariable;
        aRandomVariable= theReliabilityDomain->getRandomVariablePtr(rvNumber);
        double stdv= aRandomVariable->getStdv();
        double h= stdv/200.0;
        setCorrelationMatrix(rvNumber, 0, h);

        MatrixOperations someMatrixOperations(*correlationMatrix);

        int result= someMatrixOperations.computeCholeskyAndItsInverse();
        if (result < 0) {
                std::cerr << "XC::NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << std::endl
                        << " compute the Cholesky decomposition and its inverse " << std::endl
                        << " for the correlation matrix." << std::endl;
        }
        Matrix PerturbedInverseLowerCholesky= someMatrixOperations.getInverseLowerCholesky();
        Matrix DinverseLowerCholeskyDmean= (OrigInverseLowerCholesky - PerturbedInverseLowerCholesky) * (1/h);
        setCorrelationMatrix(0, 0, 0.0);

        // Return the final result (the four factors)
        return ( DinverseLowerCholeskyDmean * z + (*inverseLowerCholesky) * DzDmean );
}



 XC::Vector 
XC::NatafProbabilityTransformation::stdvSensitivityOf_x_to_u(Vector &x, int rvNumber)
{
        // Returns the sensitivity of 'u' with respect to [mean, stdv]
        // for the given random variable number (rvNumber)

        // Need four factors:

        // 1) (*inverseLowerCholesky)    DONE

        // 2) Vector z= x_to_z(x); 
        Vector z= x_to_z(x);

        // 3) DzDmean and DzDstdv= a vector of zeros and then:
        Vector DzDstdv(x.Size());
        NormalRV aStandardNormalRV(1,0.0,1.0,0.0); 
        RandomVariable *theRV= theReliabilityDomain->getRandomVariablePtr(rvNumber);
        if(strcmp(theRV->getType(),"NORMAL")==0)
          {
                double mu= theRV->getParameter1();
                double sigma= theRV->getParameter2();
                DzDstdv(rvNumber-1)= - (x(rvNumber-1)-mu) / (sigma*sigma);
          }
        else if(strcmp(theRV->getType(),"LOGNORMAL")==0)
          {
	    //double lambda= theRV->getParameter1();
            //double zeta= fabs(theRV->getParameter2()); // more here for negative lognormal?
                double mean= fabs(theRV->getMean()); // more here for negative lognormal?
                double stdv= theRV->getStdv();

                double a= mean*mean+stdv*stdv;
                DzDstdv(rvNumber-1)= 0.5*stdv*(log(a)-2.0*log(fabs(x(rvNumber-1))))/(pow(log(a)-2.0*log(mean),1.5)*a);
                
//                double z= ( log ( fabs(x(rvNumber-1)) ) - lambda ) / zeta; // more here for negative lognormal?
//                double e= (stdv/mean)*(stdv/mean);
//                double d= 1 + e;
//                double f= 1.0 / (mean*d*zeta);
//                DzDstdv(rvNumber-1)= stdv*((1.0-z/zeta)*f/mean);
        }
        else if (strcmp(theRV->getType(),"UNIFORM")==0) {
                double pz= 0.39894228048*exp(-z(rvNumber-1)*z(rvNumber-1)/2.0);
                double a= theRV->getParameter1();
                double b= theRV->getParameter2();
                double DzDmean= -1.0/(pz*(b-a));
                double e= -DzDmean/(b-a);
                DzDstdv(rvNumber-1)= 1.732050807*(a+b-2.0*x(rvNumber-1))*e;
        }
        else {
                std::cerr << "WARNING: Cannot compute reliability sensitivity results for " << std::endl
                        << " type of random variable number " << rvNumber << std::endl;
                DzDstdv(rvNumber-1)= 0.0; //aStandardNormalRV.getInverseCDFvalue(theRV->getCDFStdvSensitivity(x(rvNumber-1)));
        }

        // 4) The hardest part: DinverseLowerCholeskyDmean
        //                                                DinverseLowerCholeskyDstdv
        //    This is possible, but a bit tedious
        //    The procedure would be:
        //             R= L * L^T
        //            dR= dL* L^T + L * dL^T
        //            dR= dL* L^T + dL* L^T
        //            dR= 2 * dL* L^T
        //            dL= 0.5 * dR * L^(-T)
        //
        //    For now: do this by finite difference
        Matrix OrigInverseLowerCholesky= (*inverseLowerCholesky);
        RandomVariable *aRandomVariable;
        aRandomVariable= theReliabilityDomain->getRandomVariablePtr(rvNumber);
        double stdv= aRandomVariable->getStdv();
        double h= stdv/200.0;
        setCorrelationMatrix(0, rvNumber, h);

        MatrixOperations someMatrixOperations(*correlationMatrix);

        int result= someMatrixOperations.computeCholeskyAndItsInverse();
        if (result < 0) {
                std::cerr << "XC::NatafProbabilityTransformation::NatafProbabilityTransformation() - could not" << std::endl
                        << " compute the Cholesky decomposition and its inverse " << std::endl
                        << " for the correlation matrix." << std::endl;
        }
        Matrix PerturbedInverseLowerCholesky= someMatrixOperations.getInverseLowerCholesky();
        Matrix DinverseLowerCholeskyDstdv= (OrigInverseLowerCholesky - PerturbedInverseLowerCholesky) * (1/h);
        setCorrelationMatrix(0, 0, 0.0);

        // Return the final result (the four factors)
        return ( DinverseLowerCholeskyDstdv * z + (*inverseLowerCholesky) * DzDstdv );
}




int 
XC::NatafProbabilityTransformation::transform_u_to_x()
{

        Vector z= (*lowerCholesky) * (*u);
        (*x)= z_to_x(z);


        // If user has set print flag to '1' then print realization 
        if (printFlag == 1) {
                RandomVariable *theRV;
                double mean, stdv, dist; 
                char theString[80];
                sprintf(theString," CURRENT REALIZATION OF RANDOM VARIABLES:");
                std::cerr << theString << std::endl;
                for ( int i=0 ; i<nrv ; i++ )
                {
                        theRV= theReliabilityDomain->getRandomVariablePtr(i+1);
                        mean= theRV->getMean();
                        stdv= theRV->getStdv();
                        dist= ((*x)(i)-mean)/stdv; 
                        sprintf(theString," x_%d: %5.2e (%5.2f standard deviations away from the mean)",(i+1),(*x)(i),dist);
                        std::cerr << theString << std::endl;
                }
        }

        return 0;
}


int 
XC::NatafProbabilityTransformation::transform_u_to_x_andComputeJacobian()
{
        Vector z= (*lowerCholesky) * (*u);
        (*x)= z_to_x(z);


        Matrix Jzx= getJacobian_z_x((*x),z);
        (*jacobian_u_x)= (*inverseLowerCholesky) * Jzx;


        int result= theMatrixOperations->setMatrix((*jacobian_u_x));
        if (result < 0) {
                std::cerr << "XC::NatafProbabilityTransformation::transform_u_to_x() - could not set " << std::endl
                        << " the matrix in the object to perform matrix operations." << std::endl;
                return -1;
        }
        
        result= theMatrixOperations->computeInverse();
        if (result < 0) {
                std::cerr << "XC::NatafProbabilityTransformation::transform_u_to_x() - could not " << std::endl
                        << " invert Jacobian_u_x." << std::endl;
                return -1;
        }
        (*jacobian_x_u)= theMatrixOperations->getInverse();


        // If user has set print flag to '1' then print realization 
        if (printFlag == 1) {
                RandomVariable *theRV;
                double mean, stdv, dist; 
                char theString[80];
                sprintf(theString," CURRENT REALIZATION OF RANDOM VARIABLES:");
                std::cerr << theString << std::endl;
                for ( int i=0 ; i<nrv ; i++ )
                {
                        theRV= theReliabilityDomain->getRandomVariablePtr(i+1);
                        mean= theRV->getMean();
                        stdv= theRV->getStdv();
                        dist= ((*x)(i)-mean)/stdv; 
                        sprintf(theString," x_%d: %5.2e (%5.2f standard deviations away from the mean)",(i+1),(*x)(i),dist);
                        std::cerr << theString << std::endl;
                }
        }

        return 0;
}







 XC::Vector 
XC::NatafProbabilityTransformation::get_x()
{
        return (*x);
}



 XC::Vector 
XC::NatafProbabilityTransformation::get_u()
{
        return (*u);
}



 XC::Matrix 
XC::NatafProbabilityTransformation::getJacobian_x_u()
{
        return (*jacobian_x_u);
}



 XC::Matrix 
XC::NatafProbabilityTransformation::getJacobian_u_x()
{
        return (*jacobian_u_x);
}







 XC::Matrix
XC::NatafProbabilityTransformation::getJacobian_z_x(Vector x, Vector z)
{        
        RandomVariable *theRV;
        NormalRV aStandardNormalRV(1, 0.0, 1.0, 0.0);
        Matrix jacobianMatrix_z_x(nrv,nrv);
        for ( int i=0 ; i<nrv ; i++ )
        {
                theRV= theReliabilityDomain->getRandomVariablePtr(i+1);
                if (strcmp(theRV->getType(),"NORMAL")==0) {
                        double sigma= theRV->getParameter2();
                        jacobianMatrix_z_x(i,i)= 1.0 / sigma;
                }
                else if (strcmp(theRV->getType(),"LOGNORMAL")==0) {
                        double zeta= fabs(theRV->getParameter2());
                        if (x(i) == 0.0) {
                                std::cerr << "XC::NatafProbabilityTransformation::getJacobian_z_x() -- Error: value " << std::endl
                                    << "of lognormal random variable is zero in original space. " << std::endl;
                        }
                        jacobianMatrix_z_x(i,i)= 1.0 / ( zeta * x(i)  );
                }
                else {
                        double pdf= aStandardNormalRV.getPDFvalue(z(i));
                        if (pdf == 0.0) {
                                std::cerr << "ERROR: XC::NatafProbabilityTransformation::getJacobian_z_x() -- " << std::endl
                                        << " the PDF value is zero, probably due to too large step. " << std::endl;
                        }
                        jacobianMatrix_z_x(i,i)= theRV->getPDFvalue(x(i)) / pdf;
                        if (jacobianMatrix_z_x(i,i)==0.0) {
                        }
                }
        }

        return jacobianMatrix_z_x;
}




 XC::Vector
XC::NatafProbabilityTransformation::z_to_x(Vector z)
{
        RandomVariable *theRV;
        NormalRV aStandardNormalRV(1, 0.0, 1.0, 0.0);
        Vector x(nrv);
        for ( int i=0 ; i<nrv ; i++ )
        {
                theRV= theReliabilityDomain->getRandomVariablePtr(i+1);
                if (strcmp(theRV->getType(),"NORMAL")==0) {
                        double mju= theRV->getParameter1();
                        double sigma= theRV->getParameter2();
                        x(i)= z(i) * sigma + mju;
                }
                else if (strcmp(theRV->getType(),"LOGNORMAL")==0) {
                        double lambda= theRV->getParameter1();
                        double zeta= theRV->getParameter2();
                        if (zeta < 0.0) { // interpret this as negative lognormal random variable
                                x(i)= -exp ( -z(i) * zeta + lambda );
                        }
                        else {
                                x(i)= exp ( z(i) * zeta + lambda );
                        }
                }
                else {
                        x(i)= theRV->getInverseCDFvalue(  aStandardNormalRV.getCDFvalue(z(i))  );
                }
        }
        return x;
}





 XC::Vector
XC::NatafProbabilityTransformation::x_to_z(Vector x)
{
        RandomVariable *theRV;
        NormalRV aStandardNormalRV(1, 0.0, 1.0, 0.0);
        Vector z(nrv);
        for ( int i=0 ; i<nrv ; i++ )
        {
                theRV= theReliabilityDomain->getRandomVariablePtr(i+1);
                if (strcmp(theRV->getType(),"NORMAL")==0) {
                        double mju= theRV->getParameter1();
                        double sigma= theRV->getParameter2();
                        z(i)=   ( x(i) - mju ) / sigma;
                }
                else if (strcmp(theRV->getType(),"LOGNORMAL")==0) {
                        double lambda= theRV->getParameter1();
                        double zeta= theRV->getParameter2();
                        if (zeta < 0.0) { /// interpret this as a negative lognormal random variable
                                z(i)= -( log ( fabs(x(i)) ) - lambda ) / zeta;
                        }
                        else {
                                z(i)= ( log ( x(i) ) - lambda ) / zeta;
                        }
                }
                else {
                        z(i)= aStandardNormalRV.getInverseCDFvalue(theRV->getCDFvalue(x(i)));
                }
        }
        return z;
}




void XC::NatafProbabilityTransformation::setCorrelationMatrix(int pertMeanOfThisRV, int pertStdvOfThisRV, double h)
{
        // Initial declarations
        int rv1;
        int rv2;
        RandomVariable *rv1Ptr;
        RandomVariable *rv2Ptr;
        double correlation;
        double newCorrelation= 0.0;

        // Initialize correlation matrix
        correlationMatrix->Zero();

        // Put 'ones' on the diagonal
        for ( int i=0 ; i<nrv ; i++ )
        {
                (*correlationMatrix)(i,i)= 1.0;
        }

        // Get number of correlation coefficients
        int numberOfCorrelationCoefficients= 
                theReliabilityDomain->getNumberOfCorrelationCoefficients();

        // Modify each coefficient at a time and put it into the correlation matrix
        for ( int j=1 ; j<=numberOfCorrelationCoefficients ; j++ )
        {
                // Get a pointer to the correlation coefficient with tag 'j'
                CorrelationCoefficient *theCorrelationCoefficient= 
                        theReliabilityDomain->getCorrelationCoefficientPtr(j);

                // Get tags for the two involved random variables
                rv1= theCorrelationCoefficient->getRv1();
                rv2= theCorrelationCoefficient->getRv2();

                // Get value of the correlation
                correlation= theCorrelationCoefficient->getCorrelation();

                // Get pointers to the two random variables
                rv1Ptr= theReliabilityDomain->getRandomVariablePtr(rv1);
                rv2Ptr= theReliabilityDomain->getRandomVariablePtr(rv2);

                // Get the types of the two random variables
                const char *typeRv1= rv1Ptr->getType();
                const char *typeRv2= rv2Ptr->getType();

                // Compute the coefficient of variation of the random variables
                if (rv1Ptr->getMean() == 0.0  ||  rv2Ptr->getMean() == 0.0 ) {
                        std::cerr << "WARNING: The Nataf-transformation is currently implemented " << std::endl
                                << " with the assumption of non-zero mean. " << std::endl;
                }
                double cov1= rv1Ptr->getStdv() / rv1Ptr->getMean();
                double cov2= rv2Ptr->getStdv() / rv2Ptr->getMean();


                // Handle negative lognormal random variables
                // (is the closed form expressions correct for these?)
                if ( strcmp(typeRv1,"LOGNORMAL") == 0 ) {
                        if (cov1 < 0.0) {
                                cov1= cov1;
                        }
                }
                if ( strcmp(typeRv2,"LOGNORMAL") == 0 ) {
                        if (cov2 < 0.0) {
                                cov2= cov2;
                        }
                }

                // Handle mean/stdv reliability sensitivities
                // (used in finite difference schemes to obtain derivative
                //  of Cholesky decomposition)
                if (pertMeanOfThisRV !=0) {
                        if (pertMeanOfThisRV == rv1) {
                                cov1= rv1Ptr->getStdv() / (rv1Ptr->getMean()+h);
                        }
                        else if (pertMeanOfThisRV == rv2) {
                                cov2= rv2Ptr->getStdv() / (rv2Ptr->getMean()+h);
                        }
                }
                else if (pertStdvOfThisRV !=0) {
                        if (pertStdvOfThisRV == rv1) {
                                cov1= (rv1Ptr->getStdv()+h) / rv1Ptr->getMean();
                        }
                        else if (pertStdvOfThisRV == rv2) {
                                cov2= (rv2Ptr->getStdv()+h) / rv2Ptr->getMean();
                        }
                }

                // Modify the correlation coefficient according to 
                // the type of the the two involved random variables

                /////////////////////////////////////////////////////////////////////////////////
                if ( strcmp(typeRv1,"USERDEFINED") == 0  ||  strcmp(typeRv2,"USERDEFINED") == 0  ) {
                        std::cerr << " ... modifying correlation rho("<<rv1<<","<<rv2<<") for user-defined random variable..." << std::endl;
                        newCorrelation= solveForCorrelation(rv1, rv2, correlation);
                        std::cerr << " ... computed correlation for Nataf standard normal variates: " << newCorrelation << std::endl;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        double zeta2= sqrt ( log ( 1.0 + cov2 * cov2 ) );
                        newCorrelation= correlation * cov2 / zeta2;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.001 - 0.007 * cov2 + 0.118 * cov2 * cov2) * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= 1.107 * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= 1.014 * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= 1.023 * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        newCorrelation= (1.026 + 0.001*correlation - 0.178*u2
                                + 0.268*s2 - 0.001*correlation*correlation
                                + 0.178*u2*u2 - 0.679*s2*s2 - 0.003*s2*correlation)  *  correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= 1.031 * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= 1.031 * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.030 + 0.238*cov2 + 0.364*cov2*cov2) * correlation;
                }
                else if ( strcmp(typeRv1,"NORMAL") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.031 - 0.195*cov2 + 0.328*cov2*cov2) * correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        double zeta1= sqrt ( log ( 1 + pow ( cov1 , 2 ) ) );
                        newCorrelation= correlation * cov1 / zeta1;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        double zeta1= sqrt ( log ( 1.0 + cov1*cov1 ) );
                        double zeta2= sqrt ( log ( 1.0 + cov2*cov2 ) );
                        newCorrelation= 1.0 / ( zeta1 * zeta2 ) * log ( 1 + correlation * cov1 * cov2 );
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        double temp= 1.001 + 0.033*correlation + 0.004*cov1 - 0.016*cov2 
                                + 0.002*correlation*correlation + 0.223*cov1*cov1 + 0.0130*cov2*cov2;
                        newCorrelation= correlation * (temp - 0.104*correlation*cov1 + 0.029*cov1*cov2 - 0.119*correlation*cov2);
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.098 + 0.003*correlation + 0.019*cov1 + 0.025*correlation*correlation
                                + 0.303*cov1*cov1 - 0.437*correlation*cov1) * correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.011 + 0.001*correlation + 0.014*cov1 + 0.004*correlation*correlation 
                                + 0.231*cov1*cov1 - 0.130*correlation*cov1) * correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.019 + 0.014*cov1 + 0.010*correlation*correlation + 0.249*cov1*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double uu= u2*u2;
                        double ss= s2*s2;
                        double xu= correlation * u2;
                        double xs= correlation * s2;
                        double us= u2*s2;
                        double up= u2*cov1;
                        double sp= s2 * cov1;
                        double temp= 0.979 + 0.053*correlation + 0.181*u2 + 0.293*s2 + 0.002*cov1
                                - 0.004*correlation*correlation - 0.181*uu + 5.796*ss + 0.277*cov1*cov1 - 0.107*xu
                                - 0.619*xs - 0.190*correlation*cov1 - 3.976*us - 0.097*up + 0.133*sp
                                - 14.402*ss*s2 - 0.069*cov1*cov1*cov1
                                + 0.031*correlation*xs + 0.015*correlation*correlation*cov1
                                + 3.976*uu*s2 + 0.097*uu*cov1 - 0.430*ss*cov1
                                + 0.113*sp*cov1 + 1.239*correlation*us + 0.380*correlation*up;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.029 + 0.001*correlation + 0.014*cov1 + 0.004*correlation*correlation
                                + 0.233*cov1*cov1 - 0.197*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.029 - 0.001*correlation + 0.014*cov1 + 0.004*correlation*correlation
                                + 0.233*cov1*cov1 + 0.197*correlation*cov1) * correlation;
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        double temp= 1.026 + 0.082*correlation - 0.019*cov1 + 0.222*cov2 + 0.018*correlation*correlation
                                + 0.288*cov1*cov1 + 0.379*cov2*cov2;
                        newCorrelation= (temp - 0.441*correlation*cov1 + 0.126*cov1*cov2 - 0.277*correlation*cov2);
                }
                else if ( strcmp(typeRv1,"LOGNORMAL") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        double temp= 1.031+0.052*correlation+0.011*cov1-0.210*cov2+0.002*correlation*correlation+0.220*cov1*cov1+0.350*cov2*cov2;
                        newCorrelation= (temp+0.005*correlation*cov1+0.009*cov1*cov2-0.174*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= (1.001 - 0.007 * cov1 + 0.118 * cov1 * cov1) * correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        double temp= 1.001 + 0.033*correlation + 0.004*cov2 - 0.016*cov1 
                                + 0.002*correlation*correlation + 0.223*cov2*cov2 + 0.0130*cov1*cov1;
                        newCorrelation= correlation * (temp - 0.104*correlation*cov2 + 0.029*cov2*cov1 - 0.119*correlation*cov1);
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        double temp= 1.002+0.022*correlation-0.012*(cov1+cov2)+0.001*correlation*correlation+0.125*(cov1*cov1+cov2*cov2);
                        newCorrelation= (temp-0.077*(correlation*cov1+correlation*cov2)+0.014*cov1*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.104+0.003*correlation-0.008*cov1+0.014*correlation*correlation+0.173*cov1*cov1-0.296*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.014+0.001*correlation-0.007*cov1+0.002*correlation*correlation+0.126*cov1*cov1-0.090*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.023+0.000*correlation-0.007*cov1+0.002*correlation*correlation+0.127*cov1*cov1-0.000*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double uu=u2*u2;
                        double ss=s2*s2;
                        double xu=correlation*u2;
                        double xs=correlation*s2;
                        double us=u2*s2;
                        double up=u2*cov1;
                        double sp=s2*cov1;
                        double temp;
                        double xp=correlation*cov1;
                        if(correlation > 0.0) {
                                temp= 0.931+0.050*correlation+0.366*u2+0.549*s2+0.181*cov1
                                -0.055*correlation*correlation-0.515*uu+4.804*ss-0.484*cov1*cov1-0.064*xu
                                -0.637*xs+0.032*xp-4.059*us-0.319*up-0.211*sp
                                +0.052*correlation*correlation*correlation+0.227*uu*u2-10.220*ss*s2+0.559*cov1*cov1*cov1-0.042*correlation*xu
                                +0.223*correlation*xs-0.172*correlation*xp+0.028*correlation*uu+0.695*correlation*ss+0.126*correlation*cov1*cov1
                                +3.845*uu*s2+0.019*uu*cov1-1.244*us*s2+0.008*up*cov1-2.075*ss*cov1
                                +0.167*sp*cov1+0.666*correlation*us+0.386*correlation*up-0.517*correlation*sp+2.125*us*cov1;
                        }
                        else {
                                temp= 1.025+0.050*correlation-0.029*u2+0.047*s2-0.136*cov1
                                +0.069*correlation*correlation+0.178*uu+6.281*ss+0.548*cov1*cov1-0.027*xu
                                -0.686*xs+0.046*xp-3.513*us+0.231*up+0.299*sp
                                +0.063*correlation*correlation*correlation-0.226*uu*u2-17.507*ss*s2-0.366*cov1*cov1*cov1+0.051*correlation*xu
                                -0.246*correlation*xs+0.186*correlation*xp-0.001*correlation*uu+0.984*correlation*ss+0.121*correlation*cov1*cov1
                                +3.700*uu*s2+0.081*uu*cov1+1.356*us*s2+0.002*up*cov1+1.654*ss*cov1
                                -0.135*sp*cov1+0.619*correlation*us+0.410*correlation*up-0.686*correlation*sp-2.205*us*cov1;
                        }
                        newCorrelation= temp*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.031+0.001*correlation-0.007*cov1+0.003*correlation*correlation+0.131*cov1*cov1-0.132*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.031-0.001*correlation-0.007*cov1+0.003*correlation*correlation+0.131*cov1*cov1+0.132*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        double temp= 1.029+0.056*correlation-0.030*cov1+0.225*cov2+0.012*correlation*correlation+0.174*cov1*cov1+0.379*cov2*cov2;
                        newCorrelation= (temp-0.313*correlation*cov1+0.075*cov1*cov2-0.182*correlation*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"GAMMA") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        double temp= 1.032+0.034*correlation-0.007*cov1-0.202*cov2+0.000*correlation*correlation+0.121*cov1*cov1+0.339*cov2*cov2;
                        newCorrelation= (temp-0.006*correlation*cov1+0.003*cov1*cov2-0.111*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= 1.107 * correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        newCorrelation= (1.098 + 0.003*correlation + 0.019*cov2 + 0.025*correlation*correlation
                                + 0.303*cov2*cov2 - 0.437*correlation*cov2) * correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.104+0.003*correlation-0.008*cov2+0.014*correlation*correlation+0.173*cov2*cov2-0.296*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.229-0.367*correlation+0.153*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.123-0.100*correlation+0.021*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.133+0.029*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double uu=u2*u2;
                        double ss=s2*s2;
                        double xu=correlation*u2;
                        double xs=correlation*s2;
                        double us=u2*s2;
                        double temp= 1.082-0.004*correlation+0.204*u2+0.432*s2-0.001*correlation*correlation
                                - 0.204*uu+7.728*ss+0.008*xu-1.699*xs-5.338*us
                                - 19.741*ss*s2+0.135*correlation*xs+5.338*uu*s2+3.397*correlation*us;
                        newCorrelation= temp*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.142-0.154*correlation+0.031*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.142+0.154*correlation+0.031*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.109-0.152*correlation+0.361*cov2+0.130*correlation*correlation+0.455*cov2*cov2-0.728*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 )  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.147+0.145*correlation-0.271*cov2+0.010*correlation*correlation+0.459*cov2*cov2-0.467*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= 1.014 * correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        newCorrelation= (1.011 + 0.001*correlation + 0.014*cov2 + 0.004*correlation*correlation 
                                + 0.231*cov2*cov2 - 0.130*correlation*cov2) * correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.014+0.001*correlation-0.007*cov2+0.002*correlation*correlation+0.126*cov2*cov2-0.090*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.123-0.100*correlation+0.021*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.028-0.029*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.038-0.008*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double temp= 1.037-0.042*correlation-0.182*u2+0.369*s2-0.001*correlation*correlation
                                +0.182*u2*u2-1.150*s2*s2+0.084*correlation*u2;
                        newCorrelation= temp*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.046-0.045*correlation+0.006*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.046+0.045*correlation+0.006*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.036-0.038*correlation+0.266*cov2+0.028*correlation*correlation+0.383*cov2*cov2-0.229*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 )  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.047+0.042*correlation-0.212*cov2+0.353*cov2*cov2-0.136*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= 1.023 * correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        newCorrelation= (1.019 + 0.014*cov2 + 0.010*correlation*correlation + 0.249*cov2*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.023+0.000*correlation-0.007*cov2+0.002*correlation*correlation+0.127*cov2*cov2-0.000*correlation*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.133+0.029*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.038-0.008*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.047-0.047*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double temp= 1.040+0.015*correlation-0.176*u2+0.432*s2-0.008*correlation*correlation
                                +0.176*u2*u2-1.286*s2*s2-0.137*correlation*s2;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.055+0.015*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.055+0.015*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.033+0.305*cov2+0.074*correlation*correlation+0.405*cov2*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"UNIFORM") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.061-0.237*cov2-0.005*correlation*correlation+0.379*cov2*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        newCorrelation= (1.026 + 0.001*correlation - 0.178*u1
                                + 0.268*s1 - 0.001*correlation*correlation
                                + 0.178*u1*u1 - 0.679*s1*s1 - 0.003*s1*correlation)  *  correlation;

                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double uu= u1*u1;
                        double ss= s1*s1;
                        double xu= correlation * u1;
                        double xs= correlation * s1;
                        double us= u1*s1;
                        double up= u1*cov2;
                        double sp= s1 * cov2;
                        double temp= 0.979 + 0.053*correlation + 0.181*u1 + 0.293*s1 + 0.002*cov2
                                - 0.004*correlation*correlation - 0.181*uu + 5.796*ss + 0.277*cov2*cov2 - 0.107*xu
                                - 0.619*xs - 0.190*correlation*cov2 - 3.976*us - 0.097*up + 0.133*sp
                                - 14.402*ss*s1 - 0.069*cov2*cov2*cov2
                                + 0.031*correlation*xs + 0.015*correlation*correlation*cov2
                                + 3.976*uu*s1 + 0.097*uu*cov2 - 0.430*ss*cov2
                                + 0.113*sp*cov2 + 1.239*correlation*us + 0.380*correlation*up;
                        newCorrelation= temp * correlation;

                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        double ba4minus32= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus32; 
                        double s2= rv2Ptr->getStdv() / ba4minus32;
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double uu=u1*u1;
                        double ss=s1*s1;
                        double xu=correlation*u1;
                        double xs=correlation*s1;
                        double us=u1*s1;
                        double up=u1*cov2;
                        double sp=s1*cov2;
                        double temp;
                        double xp=correlation*cov2;
                        if(correlation > 0.0) {
                                temp= 0.931+0.050*correlation+0.366*u1+0.549*s1+0.181*cov2
                                -0.055*correlation*correlation-0.515*uu+4.804*ss-0.484*cov2*cov2-0.064*xu
                                -0.637*xs+0.032*xp-4.059*us-0.319*up-0.211*sp
                                +0.052*correlation*correlation*correlation+0.227*uu*u2-10.220*ss*s2+0.559*cov2*cov2*cov2-0.042*correlation*xu
                                +0.223*correlation*xs-0.172*correlation*xp+0.028*correlation*uu+0.695*correlation*ss+0.126*correlation*cov2*cov2
                                +3.845*uu*s2+0.019*uu*cov2-1.244*us*s1+0.008*up*cov2-2.075*ss*cov2
                                +0.167*sp*cov2+0.666*correlation*us+0.386*correlation*up-0.517*correlation*sp+2.125*us*cov2;
                        }
                        else {
                                temp= 1.025+0.050*correlation-0.029*u1+0.047*s1-0.136*cov2
                                +0.069*correlation*correlation+0.178*uu+6.281*ss+0.548*cov2*cov2-0.027*xu
                                -0.686*xs+0.046*xp-3.513*us+0.231*up+0.299*sp
                                +0.063*correlation*correlation*correlation-0.226*uu*u1-17.507*ss*s2-0.366*cov2*cov2*cov2+0.051*correlation*xu
                                -0.246*correlation*xs+0.186*correlation*xp-0.001*correlation*uu+0.984*correlation*ss+0.121*correlation*cov2*cov2
                                +3.700*uu*s1+0.081*uu*cov2+1.356*us*s1+0.002*up*cov2+1.654*ss*cov2
                                -0.135*sp*cov2+0.619*correlation*us+0.410*correlation*up-0.686*correlation*sp-2.205*us*cov2;
                        }
                        newCorrelation= temp*correlation;

                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double ba4minus32= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        //double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus32; 
                        double s2= rv2Ptr->getStdv() / ba4minus32;
                        double uu=u1*u1;
                        double ss=s1*s1;
                        double xu=correlation*u1;
                        double xs=correlation*s1;
                        double us=u1*s1;
                        double temp= 1.082-0.004*correlation+0.204*u1+0.432*s1-0.001*correlation*correlation
                                - 0.204*uu+7.728*ss+0.008*xu-1.699*xs-5.338*us
                                - 19.741*ss*s2+0.135*correlation*xs+5.338*uu*s1+3.397*correlation*us;
                        newCorrelation= temp*correlation;

                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double temp= 1.037-0.042*correlation-0.182*u1+0.369*s1-0.001*correlation*correlation
                                +0.182*u1*u1-1.150*s1*s1+0.084*correlation*u1;
                        newCorrelation= temp*correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double temp= 1.040+0.015*correlation-0.176*u1+0.432*s1-0.008*correlation*correlation
                                +0.176*u1*u1-1.286*s1*s1-0.137*correlation*s1;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        double ba1= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba1; 
                        double s1= rv1Ptr->getStdv() / ba1;
                        double ba2= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba2; 
                        double s2= rv2Ptr->getStdv() / ba2;
                        double o=u1;
                        double p=s1;
                        double q=u2;
                        double r=s2;
                        double u=o+q;
                        double s=p+r;
                        double oo=o*o;
                        //double pp=p*cov1;
                        //double qq=q*q;
                        double rr=r*r;
                        double oq=o*q;
                        double pr=p*r;
                        double us=oo+cov2*cov2;
                        double ss=cov1*cov1+rr;
                        double uc=oo*o+cov2*cov2*q;
                        double sc=cov1*cov1*cov1+rr*r;
                        double x=correlation;
                        double xx=correlation*correlation;
                        double temp;
                        if(correlation > 0.0) {
                                temp =1.030-0.050*x-0.056*u+0.094*s+0.009*xx-0.084*us+2.583*ss
                                        +0.100*x*u+0.010*x*s-0.485*u*s+0.270*oq-0.688*pr-0.011*xx*x
                                        +0.024*uc-10.786*sc+0.013*xx*u-0.035*xx*s+0.001*x*us-0.069*x*ss
                                        +1.174*us*s+0.004*oq*u+0.227*ss*u+2.783*pr*s+0.058*x*s*u
                                        -0.260*x*oq-0.352*x*pr-1.609*oq*s+0.194*pr*u;
                        }
                        else {
                                temp=0.939-0.023*x+0.147*u+0.668*s+0.035*xx-0.008*us+3.146*ss
                                        +0.103*x*u-0.126*x*s-1.866*u*s-0.268*oq-0.304*pr+0.011*xx*x
                                        -0.024*uc-10.836*sc-0.013*xx*u-0.035*xx*s-0.001*x*us+0.069*x*ss
                                        +1.175*us*s-0.005*oq*u-0.270*ss*u+2.781*pr*s+0.058*x*u*s
                                        -0.259*x*oq+0.352*x*pr+1.608*oq*s-0.189*pr*u;
                        }
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double temp= 1.055-0.066*correlation-0.194*u1+0.391*s1+0.003*correlation*correlation
                                +0.194*u1*u1-1.134*s1*s1+0.130*correlation*u1+0.003*correlation*s1;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double temp= 1.055+0.066*correlation-0.194*u1+0.391*s1+0.003*correlation*correlation
                                +0.194*u1*u1-1.134*s1*s1-0.130*correlation*u1-0.003*correlation*s1;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        double ba1= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba1; 
                        double s1= rv1Ptr->getStdv() / ba1;
                        double uu=u1*u1;
                        double ss=s1*s1;
                        double xu=correlation*u1;
                        double xs=correlation*s1;
                        double us=u1*s1;
                        double uq=u1*cov2;
                        double sq=s1*cov2;
                        double xq=correlation*cov2;
                        double temp=1.005 + 0.091*correlation + 0.285*u1+ 0.260*s1+ 0.199*cov2
                                - 0.023*correlation*correlation - 0.285*uu + 8.180*ss + 0.543*cov2*cov2 - 0.181*xu
                                - 1.744*xs - 0.336*xq - 5.450*us - 0.265*uq + 0.514*sq
                                -19.661*ss*s1- 0.178*cov2*cov2*cov2
                                + 0.244*correlation*xs + 0.066*correlation*correlation*cov2 - 0.001*correlation*ss
                                + 5.450*uu*s1+ 0.265*uu*cov2 - 0.986*ss*cov2
                                + 0.133*sq*cov2 + 3.488*correlation*us + 0.671*correlation*uq;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"BETA") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        const double ba4minus3= rv1Ptr->getParameter4() - rv1Ptr->getParameter3();
                        double u1= ( rv1Ptr->getMean() - rv1Ptr->getParameter3() ) / ba4minus3; 
                        double s1= rv1Ptr->getStdv() / ba4minus3;
                        double temp= 1.054+0.002*correlation-0.176*u1+0.366*s1-0.201*cov2
                                -0.002*correlation*correlation+0.176*u1*u1-1.098*s1*s1+0.340*cov2*cov2
                                -0.004*correlation*u1-0.029*s1*cov2;
                        newCorrelation= temp * correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= 1.031 * correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        newCorrelation= (1.029 + 0.001*correlation + 0.014*cov2 + 0.004*correlation*correlation
                                + 0.233*cov2*cov2 - 0.197*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.031+0.001*correlation-0.007*cov2+0.003*correlation*correlation+0.131*cov2*cov2-0.132*correlation*cov2)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.142-0.154*correlation+0.031*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.046-0.045*correlation+0.006*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.055+0.015*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double temp= 1.055-0.066*correlation-0.194*u2+0.391*s2+0.003*correlation*correlation
                                +0.194*u2*u2-1.134*s2*s2+0.130*correlation*u2+0.003*correlation*s2;
                        newCorrelation= temp * correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.064-0.069*correlation+0.005*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.064+0.069*correlation+0.005*correlation*correlation)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.056-0.060*correlation+0.263*cov2+0.020*correlation*correlation+0.383*cov2*cov2-0.332*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.064+0.065*correlation-0.210*cov2+0.003*correlation*correlation+0.356*cov2*cov2-0.211*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= 1.031 * correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        newCorrelation= (1.029 - 0.001*correlation + 0.014*cov2 + 0.004*correlation*correlation
                                + 0.233*cov2*cov2 + 0.197*correlation*cov2) * correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        newCorrelation= (1.031-0.001*correlation-0.007*cov2+0.003*correlation*correlation+0.131*cov2*cov2+0.132*correlation*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.142+0.154*correlation+0.031*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.046+0.045*correlation+0.006*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.055+0.015*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double temp= 1.055+0.066*correlation-0.194*u2+0.391*s2+0.003*correlation*correlation
                                +0.194*u2*u2-1.134*s2*s2-0.130*correlation*u2-0.003*correlation*s2;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.064+0.069*correlation+0.005*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.064-0.069*correlation+0.005*correlation*correlation)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        newCorrelation= (1.056+0.060*correlation+0.263*cov2+0.020*correlation*correlation+0.383*cov2*cov2+0.332*correlation*cov2)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE1SMALLESTVALAUE") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        newCorrelation= (1.064-0.065*correlation-0.210*cov2+0.003*correlation*correlation+0.356*cov2*cov2+0.211*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= (1.030 + 0.238*cov1 + 0.364*cov1*cov1) * correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        double temp= 1.026 + 0.082*correlation - 0.019*cov2 + 0.222*cov1 + 0.018*correlation*correlation
                                + 0.288*cov2*cov2 + 0.379*cov1*cov1;
                        newCorrelation= (temp - 0.441*correlation*cov2 + 0.126*cov2*cov1 - 0.277*correlation*cov1);                
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        double temp= 1.029+0.056*correlation-0.030*cov2+0.225*cov1+0.012*correlation*correlation+0.174*cov2*cov2+0.379*cov1*cov1;
                        newCorrelation= (temp-0.313*correlation*cov2+0.075*cov1*cov2-0.182*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.109-0.152*correlation+0.361*cov1+0.130*correlation*correlation+0.455*cov1*cov1-0.728*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.036-0.038*correlation+0.266*cov1+0.028*correlation*correlation+0.383*cov1*cov1-0.229*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.033+0.305*cov1+0.074*correlation*correlation+0.405*cov1*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        double ba2= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba2; 
                        double s2= rv2Ptr->getStdv() / ba2;
                        double uu=u2*u2;
                        double ss=s2*s2;
                        double xu=correlation*u2;
                        double xs=correlation*s2;
                        double us=u2*s2;
                        double uq=u2*cov1;
                        double sq=s2*cov1;
                        double xq=correlation*cov1;
                        double temp=1.005 + 0.091*correlation + 0.285*u2+ 0.260*s2+ 0.199*cov1
                                - 0.023*correlation*correlation - 0.285*uu + 8.180*ss + 0.543*cov1*cov1 - 0.181*xu
                                - 1.744*xs - 0.336*xq - 5.450*us - 0.265*uq + 0.514*sq
                                -19.661*ss*s2- 0.178*cov1*cov1*cov1
                                + 0.244*correlation*xs + 0.066*correlation*correlation*cov1 - 0.001*correlation*ss
                                + 5.450*uu*s2+ 0.265*uu*cov1 - 0.986*ss*cov1
                                + 0.133*sq*cov1 + 3.488*correlation*us + 0.671*correlation*uq;
                        newCorrelation= temp * correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.056-0.060*correlation+0.263*cov1+0.020*correlation*correlation+0.383*cov1*cov1-0.332*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.056+0.060*correlation+0.263*cov1+0.020*correlation*correlation+0.383*cov1*cov1+0.332*correlation*cov1)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        double rs=cov1*cov1+cov2*cov2;
                        double rc=cov1*cov1*cov1+cov2*cov2*cov2;
                        double r=cov1+cov2;
                        double pq=cov1*cov2;
                        double temp=1.086+0.054*correlation+0.104*r-0.055*correlation*correlation+0.662*rs-0.570*correlation*r+0.203*pq;
                        newCorrelation= (temp-0.020*correlation*correlation*correlation-0.218*rc-0.371*correlation*rs+0.257*correlation*correlation*r+0.141*pq*r)*correlation;
                }
                else if ( strcmp(typeRv1,"TYPE2LARGESTVALUE") == 0  &&  (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0) ) {
                        double temp= 1.065+0.146*correlation+0.241*cov1-0.259*cov2+0.013*correlation*correlation+0.372*cov1*cov1+0.435*cov2*cov2;
                        newCorrelation= (temp+0.005*correlation*cov1+0.034*cov1*cov2-0.481*correlation*cov2)*correlation;
                }
                /////////////////////////////////////////////////////////////////////////////////
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"NORMAL") == 0  ) {
                        newCorrelation= (1.031 - 0.195*cov1 + 0.328*cov1*cov1) * correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"LOGNORMAL") == 0  ) {
                        double temp= 1.031+0.052*correlation+0.011*cov2-0.210*cov1+0.002*correlation*correlation+0.220*cov2*cov2+0.350*cov1*cov1;
                        newCorrelation= (temp+0.005*correlation*cov2+0.009*cov2*cov1-0.174*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"GAMMA") == 0  ) {
                        double temp= 1.032+0.034*correlation-0.007*cov2-0.202*cov1+0.000*correlation*correlation+0.121*cov2*cov2+0.339*cov1*cov1;
                        newCorrelation= (temp-0.006*correlation*cov2+0.003*cov2*cov1-0.111*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  (strcmp(typeRv2,"EXPONENTIAL") == 0 || strcmp(typeRv2,"SHIFTEDEXPONENTIAL") == 0 ) ) {
                        newCorrelation= (1.147+0.145*correlation-0.271*cov1+0.010*correlation*correlation+0.459*cov1*cov1-0.467*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  (strcmp(typeRv2,"RAYLEIGH") == 0 || strcmp(typeRv2,"SHIFTEDRAYLEIGH") == 0 ) ) {
                        newCorrelation= (1.047+0.042*correlation-0.212*cov1+0.353*cov1*cov1-0.136*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"UNIFORM") == 0  ) {
                        newCorrelation= (1.061-0.237*cov1-0.005*correlation*correlation+0.379*cov1*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"BETA") == 0  ) {
                        const double ba4minus3= rv2Ptr->getParameter4() - rv2Ptr->getParameter3();
                        double u2= ( rv2Ptr->getMean() - rv2Ptr->getParameter3() ) / ba4minus3; 
                        double s2= rv2Ptr->getStdv() / ba4minus3;
                        double temp= 1.054+0.002*correlation-0.176*u2+0.366*s2-0.201*cov1
                                -0.002*correlation*correlation+0.176*u2*u2-1.098*s2*s2+0.340*cov1*cov1
                                -0.004*correlation*u2-0.029*s2*cov1;
                        newCorrelation= temp * correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  (strcmp(typeRv2,"TYPE1LARGESTVALUE") == 0 || strcmp(typeRv2,"GUMBEL") == 0)  ) {
                        newCorrelation= (1.064+0.065*correlation-0.210*cov1+0.003*correlation*correlation+0.356*cov1*cov1-0.211*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"TYPE1SMALLESTVALAUE") == 0  ) {
                        newCorrelation= (1.064-0.065*correlation-0.210*cov1+0.003*correlation*correlation+0.356*cov1*cov1+0.211*correlation*cov1)*correlation;
                }
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  strcmp(typeRv2,"TYPE2LARGESTVALUE") == 0  ) {
                        double temp= 1.065+0.146*correlation+0.241*cov2-0.259*cov1+0.013*correlation*correlation+0.372*cov2*cov2+0.435*cov1*cov1;
                        newCorrelation= (temp+0.005*correlation*cov2+0.034*cov1*cov2-0.481*correlation*cov1)*correlation;
                }
#ifndef _WIN32
                // VC++.net has a limit of 128 on number of nested conditionals!!!!!!!!
                else if ( (strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0 || strcmp(typeRv2,"WEIBULL") == 0)  &&  ( strcmp(typeRv2,"TYPE3SMALLESTVALUE") == 0  || strcmp(typeRv2,"WEIBULL") == 0 ) ) {
                        double temp= 1.063-0.004*correlation-0.200*(cov1+cov2)-0.001*correlation*correlation+0.337*(cov1*cov1+cov2*cov2);
                        newCorrelation= (temp+0.007*(correlation*cov1+correlation*cov2)-0.007*cov1*cov2)*correlation;
                }
#endif
                /////////////////////////////////////////////////////////////////////////////////
                else {
                        std::cerr << "Did not find the given combination of distributions in CorrelationModifier" << std::endl;
                }

                if(newCorrelation > 1.0) {
                        newCorrelation= 0.999999999;
                }
                if(newCorrelation < -1.0) {
                        newCorrelation= -0.999999999;
                }
                
                
                // Put the coefficient into the correlation matrix
                (*correlationMatrix)( ( rv1-1 ) , ( rv2-1 ) )= newCorrelation;
                (*correlationMatrix)( ( rv2-1 ) , ( rv1-1 ) )= newCorrelation;
        }

        // Here the correlation matrix should be checked for validity
        // (Whether it is close to singular or not)
        

}






double
XC::NatafProbabilityTransformation::phi2(double z_i, 
                                                                         double z_j,
                                                                         double rho)
{
        double par= z_i*z_i + z_j*z_j - 2.0*rho*z_i*z_j;

        double theExp= exp(-par/(2.0*(1.0-rho*rho)));

        double pi= 3.14159265358979;

        double result= theExp/(2.0*pi*sqrt(1.0-rho*rho));

        return result;
}



double
XC::NatafProbabilityTransformation::integrand(int rv_i,
                                                                                  double z_i, 
                                                                                  double mean_i,
                                                                                  double stdv_i, 
                                                                                  int rv_j,
                                                                                  double z_j,
                                                                                  double mean_j, 
                                                                                  double stdv_j,
                                                                                  double rho)
{
        RandomVariable *theRv_i= theReliabilityDomain->getRandomVariablePtr(rv_i);
        RandomVariable *theRv_j= theReliabilityDomain->getRandomVariablePtr(rv_j);
        NormalRV aStandardNormalRV(1,0.0,1.0,0.0); 

        double x_i= theRv_i->getInverseCDFvalue(aStandardNormalRV.getCDFvalue(z_i));
        double x_j= theRv_j->getInverseCDFvalue(aStandardNormalRV.getCDFvalue(z_j));

        double thePhi2= phi2(z_i, z_j, rho);

        return ( (x_i-mean_i)/stdv_i 
                   * (x_j-mean_j)/stdv_j 
                   * thePhi2 );
}


double
XC::NatafProbabilityTransformation::doubleIntegral(int rv_i,
                                                                                  double mean_i,
                                                                                  double stdv_i, 
                                                                                  int rv_j,
                                                                                  double mean_j, 
                                                                                  double stdv_j,
                                                                                  double rho)
{
        // The grid of integration points:
        // 1, 2, ..., i, ..., 2*n  in x-direction with intervals h
        // 1, 2, ..., j, ..., 2*m  in y-direction with intervals k
        int i, j;
        

        // Selected integration boundaries
        double z_i0= -5.0;
        double z_in=  5.0;
        double z_j0= -5.0;
        double z_jm=  5.0;

        // Half the number of integration points
        int n= 100;
        int m= 100;
        
        // Integration point intervals
        double h= (z_in-z_i0)/(2.0*n);
        double k= (z_jm-z_j0)/(2.0*m);

        // Grid of integration points
        Vector z_i(2.0*n);
        Vector z_j(2.0*m);
        for (i=1; i<=2*n; i++) {
                z_i(i-1)= z_i0 + (i-1)*h;
        }
        for (j=1; j<=2*m; j++) {
                z_j(j-1)= z_j0 + (j-1)*k;
        }

        // Computing sums (naming terms according to p. 126 in "Numerical Methods" by Faires & Burden)
        double term1= 0.0;
        double term2= 0.0;
        double term3= 0.0;
        double term4= 0.0;
        double term5= 0.0;
        double term6= 0.0;
        double term7= 0.0;
        double term8= 0.0;
        double term9= 0.0;
        double term10= 0.0;
        double term11= 0.0;
        double term12= 0.0;
        double term13= 0.0;
        double term14= 0.0;
        double term15= 0.0;
        double term16= 0.0;

        // No sum terms
        term1= integrand( rv_i, z_i( 0        ) , mean_i, stdv_i, 
                                   rv_j, z_j( 0        ) , mean_j, stdv_j, rho);
        term4= integrand( rv_i, z_i( 2*n    -1) , mean_i, stdv_i, 
                                   rv_j, z_j( 0        ) , mean_j, stdv_j, rho);
        term13= integrand(rv_i, z_i( 0        ) , mean_i, stdv_i, 
                                   rv_j, z_j( 2*m    -1) , mean_j, stdv_j, rho);
        term16= integrand(rv_i, z_i( 2*n    -1) , mean_i, stdv_i, 
                                   rv_j, z_j( 2*m    -1) , mean_j, stdv_j, rho);

        // Single sum over n terms
        for (i=1; i<=n; i++) {
                term2 += integrand(rv_i, z_i( 2*i    -1) , mean_i, stdv_i, 
                                       rv_j, z_j( 0        ) , mean_j, stdv_j, rho);
                term3 += integrand(rv_i, z_i( 2*i-1  -1) , mean_i, stdv_i, 
                                       rv_j, z_j( 0        ) , mean_j, stdv_j, rho);
                term14 += integrand(rv_i, z_i( 2*i    -1) , mean_i, stdv_i, 
                                        rv_j, z_j( 2*m    -1) , mean_j, stdv_j, rho);
                term15 += integrand(rv_i, z_i( 2*i-1    -1) , mean_i, stdv_i, 
                                        rv_j, z_j( 2*m      -1) , mean_j, stdv_j, rho);
        }
        term2 -= integrand(rv_i, z_i( 2*n  -1) , mean_i, stdv_i, 
                                   rv_j, z_j( 0      ) , mean_j, stdv_j, rho);
        term14 -= integrand(rv_i, z_i( 2*n    -1) , mean_i, stdv_i, 
                                rv_j, z_j( 2*m    -1) , mean_j, stdv_j, rho);

        // Single sum over m terms
        for (j=1; j<=m; j++) {
                term5 += integrand(rv_i, z_i( 0        ) , mean_i, stdv_i, 
                                       rv_j, z_j( 2*j    -1) , mean_j, stdv_j, rho);
                term8 += integrand(rv_i, z_i( 2*n    -1) , mean_i, stdv_i, 
                                       rv_j, z_j( 2*j    -1) , mean_j, stdv_j, rho);
                term9 += integrand(rv_i, z_i( 0          ) , mean_i, stdv_i, 
                                       rv_j, z_j( 2*j-1    -1) , mean_j, stdv_j, rho);
                term12 += integrand(rv_i, z_i( 2*n     -1) , mean_i, stdv_i, 
                                        rv_j, z_j( 2*j-1   -1) , mean_j, stdv_j, rho);
        }
        term8 -= integrand(rv_i, z_i( 2*n    -1) , mean_i, stdv_i, 
                               rv_j, z_j( 2*m    -1) , mean_j, stdv_j, rho);
        
        // Double sum terms
        for (j=1; j<=(m-1); j++) {
                for (i=1; i<=(n-1); i++) {
                        term6 += integrand(rv_i, z_i( 2*i    -1) , mean_i, stdv_i, 
                                                           rv_j, z_j( 2*j    -1) , mean_j, stdv_j, rho);
                }
        }
        for (j=1; j<=(m-1); j++) {
                for (i=1; i<=(n); i++) {
                        term7 += integrand(rv_i, z_i( 2*i-1    -1) , mean_i, stdv_i, 
                                                           rv_j, z_j( 2*j      -1) , mean_j, stdv_j, rho);
                }
        }
        for (j=1; j<=(m); j++) {
                for (i=1; i<=(n-1); i++) {
                        term10 += integrand(rv_i, z_i( 2*i       -1) , mean_i, stdv_i, 
                                                            rv_j, z_j( 2*j-1     -1) , mean_j, stdv_j, rho);
                }
        }
        for (j=1; j<=(m); j++) {
                for (i=1; i<=(n); i++) {
                        term11 += integrand(rv_i, z_i( 2*i-1     -1) , mean_i, stdv_i, 
                                                            rv_j, z_j( 2*j-1     -1) , mean_j, stdv_j, rho);
                }
        }


        double par1= term1 + 2.0*term2 + 4.0*term3 + term4;
        
        double par2= term5 + 2.0*term6 + 4.0*term7 + term8;

        double par3= term9 + 2.0*term10 + 4.0*term11 + term12;

        double par4= term13 + 2.0*term14 + 4.0*term15 + term16;

        double result= h*k/9.0 * (par1 + 2.0*par2 + 4.0*par3 + par4);

        return result;
}



double
XC::NatafProbabilityTransformation::residualFunction(double rho_original, 
                                                                                                 double rho,
                                                                                                 double rv_i, 
                                                                                                 double mean_i, 
                                                                                                 double stdv_i, 
                                                                                                 double rv_j, 
                                                                                                 double mean_j, 
                                                                                                 double stdv_j)
{
        double result= rho_original - doubleIntegral(rv_i, mean_i, stdv_i, rv_j, mean_j, stdv_j, rho);

        return result;
}




double
XC::NatafProbabilityTransformation::solveForCorrelation(int rv_i, int rv_j, double rho_original)
{
        RandomVariable *theRv_i= theReliabilityDomain->getRandomVariablePtr(rv_i);
        RandomVariable *theRv_j= theReliabilityDomain->getRandomVariablePtr(rv_j);

        double mean_i= theRv_i->getMean();
        double mean_j= theRv_j->getMean();

        double stdv_i= theRv_i->getStdv();
        double stdv_j= theRv_j->getStdv();

        double result= 0.0;

        double tol= 1.0e-6;
        double pert= 1.0e-4;

        double rho_old= rho_original;
        double rho_new;
        double f;
        double df;
        double perturbed_f;

        for (int i=1;  i<=100;  i++ )  {

                // Evaluate function
                f= residualFunction(rho_original,
                                         rho_old,
                                                         rv_i, mean_i, stdv_i, 
                                                         rv_j, mean_j, stdv_j);

                // Evaluate perturbed function
                perturbed_f= residualFunction(rho_original,
                                                                           (rho_old+pert),
                                                                           rv_i, mean_i, stdv_i, 
                                                                           rv_j, mean_j, stdv_j);

                // Evaluate derivative of function
                df= ( perturbed_f - f ) / pert;

                if ( fabs(df) < 1.0e-15) {
                        std::cerr << "WARNING: XC::NatafProbabilityTransformation::solveForCorrelation() -- " << std::endl
                                << " zero derivative in Newton algorithm. " << std::endl;
                }
                else {

                        // Take a Newton step
                        rho_new= rho_old - f/df;
                        
                        // Check convergence; quit or continue
                        if (fabs(1.0-fabs(rho_old/rho_new)) < tol) {
                                result= rho_new;
                                return result;
                        }
                        else {
                                if (i==100) {
                                        std::cerr << "WARNING: XC::NatafProbabilityTransformation::solveForCorrelation() -- " << std::endl
                                                << " Newton scheme did not converge. " << std::endl;
                                        result= 0.0;
                                        return result;
                                }
                                else {
                                        rho_old= rho_new;
                                }
                        
                        }
                }
        }

        return result;

}













